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Abstract. An eigcnfunction expansion method involving hypergeometric functions is used to solve the 
partial differential equation governing the transport of radiation in an X-ray pulsar accretion column 
containing a radiative shock. The procedure yields the exact solution for the Green's function, which 
describes the scattering of monochromatic radiation injected into the column from a source located near the 
surface of the star. Collisions between the injected photons and the infalling electrons cause the radiation 
to gain energy as it diffuses through the gas and gradually escapes by passing through the walls of the 
column. The presence of the shock enhances the energization of the radiation and creates a power-law 
spectrum at high energies, which is typical for a Fermi process. The analytical solution for the Green's 
function provides important physical insight into the spectral formation process in X-ray pulsars, and it 
also has direct relevance for the interpretation of spectral data for these sources. Additional interesting 
mathematical aspects of the problem include the establishment of a closed-form expression for the quadratic 
normalization integrals of the orthogonal eigenfunctions, and the derivation of a new summation formula 
involving products of hypergeometric functions. By taking various limits of the general expressions, we also 
develop new linear and bilinear generating functions for the Jacobi polynomials. 



I. INTRODUCTION 

In this article, methods of classical analysis are employed to obtain the exact solution for the Green's 
function describing the spectrum of radiation emitted by an X-ray pulsar. Beyond the direct physical 
relevance of the Green's function, the method of solution also yields several additional results of mathematical 
interest, including a new summation formula involving products of two hypergeometric functions, as well as 
new linear and bilinear generating functions for the Jacobi polynomials. We also obtain an exact expression 
for the quadratic normalization integrals of the orthogonal hypergeometric eigenfunctions. Before proceeding 
with the main derivation, some physical background is called for. The radiation produced in bright X-ray 
pulsars is powered by the gravitational accretion (inflow) of ionized gas that is channeled onto the poles of a 
rotating neutron star by the strong magnetic field. In these sources, the radiation pressure greatly exceeds 
the gas pressure, and therefore the pressure of the photons governs the dynamical structure of the accretion 
flow. It follows that the gas must pass through a radiation-dominated shock on its way to the stellar surface, 
and the kinetic energy of the gas is carried away by the high-energy radiation that escapes from the column!^ 
The strong gradient of the radiation pressure decelerates the material to rest at the surface of the star, and 
the compression of the infalling gas drives its temperatures up to a few million Kelvins. The gas therefore 
radiates X-rays, which appear to pulsate due to the star's spin. However, the observed X-ray spectrum is 
nonthermal, indicating that nonequilibrium processes are playing an important role in the formation of the 
radiation distribution. 

The nonthermal shape of the spectrum is primarily due to the flow compression, which causes Fermi 
energization of the photons as they collide with infalling electrons in the column, until the radiation escapes 
from the column into space. Our primary goal in this article is to obtain an exact solution for the Green's 
function describing the upscattering of soft, monoenergetic photons injected by a source located in the base of 
the accretion column, near the surface of the star. The Green's function contains a complete representation 
of the fundamental physics governing the propagation of the photons in the physical and energy spaces. 
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Since the transport equation governing the radiation distribution is linear, we can compute the solution 
associated with an arbitrary source distribution via convolution. Hence the Green's function provides the 
most direct means for exploring the relationship between the physics occurring in the accretion shock and 
the production of the observed nonthermal X-radiation. 

II. FUNDAMENTAL EQUATIONS 

We assume that the accretion column is cylindrical, and we define x as the spatial coordinate measured 
along the column axis. The gas flows through the column onto the stellar surface with velocity v. We define 
the Green's function, f G (cc , x, e , e), as the radiation distribution at location x and energy e resulting from 
the injection of Nq photons per second with energy eo from a monochromatic source at location xq inside 
the column. In a steady-state situation, f G satisfies the transport equatior0El 

v df G = & £ ££/g + _££ ( c d/ G \ + N S{e- e Q )S(x ~ x ) f a ^..j/ % j ,y 

dx dx 3 de dx \3n e cr^ dx J 7r7 "o e o ^ csc 

where n e is the electron number density, crn is the electron scattering cross section for photons propagating 
parallel to the x-axis, tq is the radius of the column, vq is the flow speed at the source location, c is the 
speed of light, and t csc is the mean time photons spend in the column before escaping through the walls into 
space. The total radiation number and energy densities associated with the distribution function f G are, 
respectively, 

/>oo />oo 

n G (x)= e 2 f G de, U Q (x) = / e 3 f G de. (2) 
Jo Jo 

The terms in Q represent, from left to right, the comoving (convective) time derivative, first-order Fermi 

energization ("bulk Comptonization" ) of the radiation in the converging flow, spatial diffusion of the photons 

parallel to the column axis, the monochromatic photon source, escape of radiation from the column, and the 

possible absorption of radiation at the source location, respectively. In physical terms, the first-order Fermi 

energization corresponds to the PdV work done on the radiation by the compression of the background 

plasma as it accretes onto the stellar surfacePI The dimensionless constant (3 expresses the strength of the 

absorption (if any) occurring at the source location, and the mean escape time is given by 

r§ n e <r± . . 

^esc = j \y) 

C 

where a± is the electron scattering cross section for photons propagating perpendicular to the column axis. 
In general, crn ^ <jj_ due to the influence of the strong magnetic field, which is directed parallel to the axis 
of the column. Absorption at the source location is expected if the photons are produced in a blackbody 
"mound" of dense gas near the base of the accretion columnP because a perfect blackbody acts as both a 
source and a sink of radiationEl 

The flux of electrons flowing down the column is denoted by J = n e v. In our cylindrical, steady-state 
problem, J maintains a constant value. BeckeJ^ demonstrated that in order for the inflowing matter to 
come to rest at the stellar surface as required, the parameters tq, J, eni, and a± must satisfy the dynamical 
constraint 

r 2 J 2 a ± a ]l = lc 2 . (4) 

In general, radiation-dominated shocks are continuous velocity transitions, with an overall thickness of a few 
Thomson scattering lengths, unlike standard (discontinuous) gas-mediated shocks.^ The exact solution for 
the inflow velocity v as a function of the spatial coordinate x is given bj^^ 
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where v c is the flow velocity at the sonicpoint, which is related to the stellar mass M*, the stellar radius 
and the gravitational constant G via^l 
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The quantity x s t appearing in (J5J is the distance between the sonic point and the stellar surface, which can 
be evaluated using Eq. (4.16) from Ref. 4 to obtain 

/ \ 1/2 



**- ma) -u;- (7) 

According to JSJl, the flow does come to rest at the surface of the star as required, since v(x st ) = 0. 
Furthermore, the constancy of the electron flux J in our cylindrical, steady-state problem implies that the 
electron number density n e is a function of x because v varies with the height inside the column [see Eq.lJSj]. 
Further simplification is possible if we work in terms of the new spatial variable y, defined by 

y(x) = f ~ J . (8) 

Note that y — > in the far upstream region (x — > — oo), and y — > 1 at the surface of the star (x — > x st ). 
Based on ^ and © , we find that the variation of the velocity v as a function of the new variable y is given 
by the simple expression 

v -^ = W-y)- (9) 

v c 4 

By combining Q, and © with the derivative relation 

1/2 

y- 1 , (io) 



dy 2^3 V^ll 

we can transform the transport equation Q for f G from x to y to obtain 



3Pv S(y - y Q ) f G _ 3 N S(e - e ) S(y - y ) 
7v c 7 it rl v c 



(11) 



where yo = y(xo) denotes the value of y at the source location. According to ©, the flow velocity at the 
source, Vo, is related to v c and yo by 

Note that we can write the Green's function as either f G (xo,x, eo,e) or f G (i/o,y,£o,e) since the parameters 
(x,xq) and (y,yo) are interchangeable via JSJ. 

III. SOLUTION FOR THE GREEN'S FUNCTION 

The physical model considered here includes Fermi energization, which tends to boost the energy of the 
injected photons as they collide with high-energy electrons streaming down through the accretion column 
towards the surface of the neutron star. Moreover, since no process that can lower the photon energy is 
included in the model, all of the photons injected from a source of monochromatic radiation with energy 
e = eo must at later times have energy e > cq. It follows that f G = for e < cq. When e > eo, (|11|) is 
separable in energy and space using the functions 

h(e,y) = e- x g(X,y) , (13) 
where A is the separation constant, and the spatial function g satisfies the differential equation 
/, ,d 2 g (l-hy\dg (\y + y-l\ 3f3v 5(y-y ) 

y{1 - V) d^ + {—) d~y + {^y—) 9 TV C 9 ■ (U) 

In order to avoid an infinite spatial diffusion flux at y = yo, the function g must be continuous there, and 
consequently we obtain the condition 



A[ff(A,y)] 



lirna(A,y + e) - g{\vo - e) = 0. (15) 

e— >0 



We can also derive a jump condition for the derivative dg/dy at the source location by integrating i|14fl with 
respect to y in a small region around y — yo- The result obtained is 



dy 



|£ 5 (A,y„), (16) 



where we have used <|12ll to substitute for v . 

The homogeneous version of l|14|) obtained when y ^ yo has fundamental solutions given by 

¥>i(Ky) = yF{a, b; c;y) , (17) 

<pl(\,y)=y- 1 / 4 F(a-5/4, 6 - 5/4 ; 2 - c ; y) , (18) 
where F(a, 6; c; z) denotes the hypergeometric functionpl and the parameters a, 6, and c are defined by 



9- V17+16A . 9 + V17+16A 9 
a = , o = , c = — , 19 
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and therefore a + b = 



3.1 Asymptotic Analysis 

The source photons injected into the flow are unable to diffuse very far upstream due to the high speed of 
the inflowing electrons. Most of the photons escape through the walls of the column within a few scattering 
lengths of the source, and therefore we conclude that the function g must vanish in the upstream limit. 
y — > 0. Asymptotic analysis indicates that the function ipi (A, y) — > in the limit y — > as required, but 
ip*(X,y) diverges and therefore it cannot be utilized in the upstream region (y < yo). Hence g must be 
given by tp± for y < yo- Conversely, in the downstream limit, the gas settles onto the surface of the star and 
therefore g should approach a constant as y — > 1. These conditions are satisfied if A is equal to one of the 
eigenvalues, A„, which are associated with the spatial eigenfunctions, g n (y), defined by 

g n (y) =g(X n ,y) . (20) 

In order to obtain a complete understanding of the global behavior of the eigenfunctions, we must also 
consider the asymptotic behaviors of the two functions (fx and (f\ m the downstream region, which are 
discussed below. 

The hypergeometric functions appearing in l|17|) and (|18|) can be evaluated at y = 1 using Eq. (15.1.20) 
from Abramowitz & Stegun^ which gives for general values of a, 6, and c 

„/ , n r(c) Tie — a — b) 

F(a,b;c;l = ; \ -( . 21 

1 (c — a) I (c — b) 

However, for the values of a, b, and c in (|17|) and l|18|) . we find that [see Eq. (|19fl ] 

c- o-6 = 0, (22) 

and therefore the hypergeometric functions F(a , 6 ; c; y) and Fia — 5/4 , b — 5/4 ; 2 — c ; y) each diverge in 
the downstream limit y — ► 1. Since the eigenfunction g n should approach a constant as y — > 1 based on 
physical considerations, we conclude that in the downstream region (y > yo), g n must be represented by a 
suitable linear combination of ipi and tp\ that remains finite as y — ► 1. In order to make further progress, 
we need to employ Eq. (15.3.10) from Abramowitz & Stegunpl which yields for general a, b, and y 



T-i/ j | j \ T(a + 6) ^ ( a )„(6)„ 
F(a ) 6;a + 6;y) = — — - ^ 



2<f(n+ 1) - *(a + n) 



+ n) - ln(l - y) 



(1 - y)" , (23) 



where 



(24) 



Asymptotic analysis of this expression reveals that in the limit y — ► 1, the logarithmic divergences of the 
two functions ipi and p\ can be balanced by creating the new function 

MW) - T^f^ V*\V) - r( ^-; } (o) ri(A,y) , (25) 

which remains finite as y — > 1. Hence <£>2 represents the fundamental solution for g n in the region downstream 
from the source. We can use the asymptotic behaviors of ipi and (p* to show that 

. 7T fcot(7ra) + cot(7r6)l ,„„. 
lim ^(A,y) = 1 r ' 1 , n ■ 26 
T(aj r(l - 6) 

Since the solutions and ^2 are applicable in the upstream and downstream regions, respectively, the 
global expression for the eigenfunction g n is therefore given by 

9n(y) = \ (27) 

[-B n y> 2 (A n! y) , y>y , 

where the constant B n is evaluated using the continuity condition [Eq. I|15f) ] , which yields 

It follows from l|26l) , 127(1 , and l|28() that the downstream value of g n is given by 

Um = 7r[cot(7r ) +cot(7r&)] gr^p) 

y-i r(a)r(l-6) <p 2 (A n ,y ) 

Conversely, in the upstream region, ipi — > y, and therefore we have the asymptotic behavior 

U m ^iM = 1 . (30) 

w->0 y 



3.2 Eigenvalue Equation 



We can combine ((16f) . I|27|l. and (|28jl to show that the eigenvalues A n satisfy the equation 

M/ ,, \ 3/3^i(A Tl ,y )y)2(A ra ,y ) 

VK(A„,y ) : =0, (31) 

4 2/0 

where the Wronskian of the two functions <f\ and <fi is defined for general values of A and y by 

Further progress can be made by deriving an analytical expression for the Wronskian. We begin by writing 
the differential equation l|14[) governing the two functions (pi and ip2 in the self-adjoint form 



d 
dy 



1/4/1 \ dL P 



X 



4„3/4 



ip — T ip = , 



where 



T 



l-y 3Pv 6(y- y ) 



(33) 



(34) 



4y 7 /4 7v c y 3 / 4 

By applying (l-i-il) to the function <f2 and multiplying the result by ip±, and then subtracting from this the 
same expression with <fi and if2 interchanged, we obtain 



which can be rewritten as 



dy 



1/4/, \ d(f2 

y 1 0--y)-r- 

dy 



f>2 



dy 



1/4 f-, \ "V 1 

y 1 (1-2/) 



dip 1 
dy 







y V4 (1 _ y) ^ + ^ [ y l/4 (1 _ y) 



= 



(35) 
(36) 



where we have made use of the result 



— = d2{p2 - <p 2 (37) 
dy dy 2 dy 2 



Equation can rearranged in the form 

d\nW _ d_ 
dy dy 

which can be integrated to obtain the exact solution 



1/4 (1-1/) , (38) 



where -D(A) is an integration constant that depends on A but not on y. The exact dependence of D on A 
can be derived by analyzing the behaviors of the functions tpi and (p2 in the limit y — > 0. For small values 
of y, we have the asymptotic expressions^ 

¥>i -* V > y -> , 

Combining (|32[) and l|4U|) . we find that asymptotically, 

^IrTW^ 1/4 ' ^ 

4 I (a) 1 (2 — c) 
Comparing this result with l|39fl . we conclude that 

and therefore the exact solution for the Wronskian for general values of A and y is given by 

WM = l Igzf) f^. (43, 
4 r(a)r(2 - c) 1 - y 

Substituting for W in l|31(l using (|43[1 . we can rewrite the eigenvalue equation in the equivalent form 

5 r(l — a) y 3 ^ 

o r/^r/n nT" 2 =P<Pi(*n,yo)<P2(*n,yo) , (44) 

3 F(o) T(2 - c) 1 - yo 

where a and 6 are functions of A„ by virtue of (|19|1 . and c = 9/4. The roots of this expression are the 
eigenvalues A„, and the associated eigenfunctions are evaluated using l(2T|) . The first eigenvalue, Ao, is 
especially important because it determines the power-law shape of the high-energy portion of the Green's 
function [see Eq. l|T3)l]. 

In Figure 1 we plot the first eigenvalue Ao as a function of the dimcnsionless parameters (3 and yo- Note 
that Ao is a double-valued function of yo for fixed f3, which is a consequence of the imposed velocity profile 
[Eq. 10]. Physically, this behavior reflects the fact that it is always possible to achieve a desired amount 
of compression (first-order Fermi energization) by placing the source in a specific location in either the 
upstream or downstream regions of the flow. We also observe that if we increase the absorption parameter 
f3 while holding j/o fixed, then Ao increases monotonically, and therefore the high-energy spectrum becomes 
progressively steeper. This behavior is expected physically because as the absorption parameter is increased, 
the injected photons spend less time on average being energized by collisions with electrons before either 
escaping from the column or being absorbed at the source location. The decreased amount of energization 
naturally leads to a steepening of the radiation spectrum. When f3 = 0, no absorption occurs, and the index 
Ao achieves it minimum (limiting) value of 4. This limit is, however, unphysical since it yields a divergent 
result for the total photon energy density U G according to J5J). Nonetheless, the case with (3 = is interesting 
from a mathematical viewpoint, and for that reason it is further discussed in section V. 
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Figure 1. First eigenvalue Ao of the Green's function expansion plotted as a function of 
the source location yo for the indicated values of the absorption parameter j3. Note the 
steepening of the radiation spectrum that occurs when (3 is increased for a fixed value of yo, 
which reflects the decreasing residence time for the photons in the plasma (see the discussion 
in the text). 



3.3 Orthogonality of the Eigenfunctions 



We shall next demonstrate that the eigenfunctions g n (y) form an orthogonal set, which is an extremely 
useful property. This is a standard Sturm-Liouville problem and therefore we follow the usual procedure. 
Let us suppose that g n (y) and g m (y) are two eigenfunctions corresponding to the distinct eigenvalues X n 
and A m , respectively The functions g n and g m each satisfy the differential equation 1)140. and therefore we 
can utilize the self-adoint form to write [cf. Eq. (|33|) ] 



g« 



d 

dy 



y 1/4 (i-y) 



dgn 
dy 



A, 



and 



./ 4(1 _„,^ 



+ 



4 y 3 /4 

A m 

4u 3 / 4 



g n - Tg r , 



9r> 



Tg n 



, 



, 



(45) 



(46) 



where T is given by l|34l) . Subtracting the second equation from the first yields, after integrating by parts 
with respect to y from y = to y = 1, 

l 

dg m dg n 



(Xn-Xm) / y~ 3/4 <?„ (tf) g m (y) dy = 4 y 1 / 4 (1 - y) 



gn 



dy 



dy 



(47) 



Based on the asymptotic behaviors of the eigenfunctions g n and g m given by i|29|) and i|30|l , we find that the 
right-hand side of l|47J) vanishes exactly, and therefore we obtain 

(A„ - A m ) f y- 3 ^ g n (y) g m (y) dy = , (48) 
Jo 

which establishes the orthogonality of the eigenfunctions. The set of eigenfunctions is also complete according 
to the Sturm-Liouville theorem. Since the eigenfunctions are orthogonal, the Green's function can be 
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expressed as the infinite series 

oo / \ — A„ 

fa (2/0, V, eo, e) = V A n I — ) 5 „(y) , (49) 

for e > eo, where the expansion coefficients A n are computed by employing the orthogonality of the 
eigenfunctions along with the condition 

12 N 



/o(j/o,J/,eo,e) 



- 2/o) , (50) 



which is obtained by integrating the transport equation l|lf|) with respect to e in a small range surrounding 
the injection energy eo- The result obtained for the nth expansion coefficient is 

1 2 N y- 3/i g n (y ) 

7 * r £ V c £n 

where the quadratic normalization integrals, £„, are defined by 



An = -;_»^o (gl) 



[\r 3/4 g 2 n (y)dy 

Jo 



(52) 



As an alternative to numerical integration, in section 3.4 we derive a closed-form expression for evaluating 
the normalization integrals based directly on the associated differential equation. 



3.4 Quadratic Normalization Integrals 



The direct computation of the normalization integrals £„ via numerical integration is costly and time 
consuming, and therefore it is desirable to have an alternative procedure available for their evaluation. In 
fact, it is possible to derive an analytical expression for the normalization integrals based on manipulation 
of the fundamental differential equation (|14|l governing the eigenfunctions g n (y). 

Let us suppose that g(X,y) is a general solution to (|14|l for an arbitrary value of A (i.e., not necessarily 
an eigenvalue) with the asymptotic (upstream) behavior 

9{\y)->y> y^o, (53) 

which is the same as the upstream behavior of the eigenfunction g n (y) [see Eq. H30fl ]. We also stipulate that 
g must be continuous at y = j/o; an d that it satisfies the derivative jump condition given by (I16|) . After a 
bit of algebra, we find that the global solution for g consistent with these requirements can be expressed as 



g{\y) = 



</>i(A,y) 



y < yo 



Jl + 6)ipi(A,y) + b(p 2 (X,y) , y > y , 

where the coefficients a and b are given by 

- _ 3f3<Pi(\yo)<P2{\yo) f 3/3(pl(\,yo) 

' Ay Q W(X,y ) ' 



(54) 



(55) 



Ay Q W(X,y ) 

and the Wronskian W is evaluated using (|43(l . 

Comparing the general solution for g(X,y) with the solution for the eigenfunction g n {y) given by l|27|) . we 
note that 

(56) 



lim a = — 1 



lim b = B n 

A^A„ 



We can now use the self-adoint form of ljl4|l to write [cf. Eqs. Q45J) and (|46Jl ] 



0_ 

Oy 



/ 1/4 (i - y) 



dg_ 

dy 



X 



and 



4 y 3 / 4 
X n 



g-Tg}=0. 



t n s, — n 



(57) 



where T is defined by l|34fl . Subtracting the second equation from the first and integrating by parts from 
y = to y = 1 yields 

1 

dg n , v dg 



(A - A n ) / y- 3 / 4 g(\,y)g n (y)dy = 4y 1 ^(l-y) 



g(\ y)-; — gn{y) 5- 



(59) 



Since g — ► y and g„ — > y as y — > 0, we conclude that the evaluation at the lower bound y — on the 
right-hand side yields zero, and consequently in the limit A — * A„ we obtain for the quadratic normalization 
integral € n [see Eq. 



IT — 



L - 3 /4 2 ( w r V /4 (1 - y) [g(A, y) (<W^) - g n (y) (dg/dy)} 







A^A„ 



A — A n 



(60) 



j/=i 



The numerator and denominator on the right-hand side of (|60J) each vanish in the limit A — * X n , and 
therefore we can employ L'Hopital's rule to show that (e.g., BeckcJ^f) 



<£„= lim 4y 1 / 4 (l-y) 



A— >A 



dg_dgn_ d 2 g 
dy dy " <9ydA 



(61) 



Substituting the analytical forms for g n (y) and y(A, y) given by H27[) and (|54|l . respectively, we find that (1611) 
can be rewritten as 



£n= lim 4y 1 / 4 (l-y)E„ 

1 



^ (A ' ~d~X ~~dX~ "dy ^ 2(A ' y) 3^A 



(62) 



A=A„ 



where we have also utilized (|32ll and Q56J1. Based on the asymptotic behavior of (p% [see (I26ll ]. we conclude 
that the final two terms on the right-hand side of Ij62|) contribute nothing in the limit y — * 1, and therefore 
our expression for € n reduces to 



£„= lim 4y 1 ^(l~y)B n W(X,y)^ 



(63) 



A=A„ 



Since y = 1 is a singular point of the differential equation (|14l) , it is convenient to employ the relation [see 
Eq. 03 

W(X, y) y 1 / 4 (1 - y) = W(X, y ) y^ 4 (1 - y ) , (64) 
which allows us to transform the evaluation in 163|) from y = 1 to y = yo to obtain the equivalent result 



<f a 1 / 4 fi \~v<rr\ ,fi(X n ,y )d\na 

£n = 4y / {l-y )aW{X,y ) — tt r — 77- 

<MA„,y ) dX 



(65) 



A=A„ 



where we have also substituted for B n using 
using H55f) . which yields 

din a dlnifi dhnp2 dhxW 



dX dX dX dX ' 

where the derivative of the Wronskian is given by [see Eqs. I|19fl and (|43[) 

dh\W *(a) + *(l-a) 



The derivative on the right-hand side can be evaluated 

(66) 



and 



dX 
9(z) 



(i7 + i6A)V2 ' 
i dr(z) 



(67) 



(68) 



T(z) dz ■ 

Combining (JSHJl, and ljfi7|) . we find that that the quadratic normalization integrals can be evaluated 

using the closed-form expression 

€ n = K(X n ,y ) , (69) 

where 

'*(a) + *(! - a) dlnipi dlntp 2 ' 



K(X,y)^3f3y~V 4 (l-y)ri(X,y) 



(17+16A) 1 /2 



OX 



OX 



(70) 




Figure 2. Green's function f a (yo,y,eo,e) [Eq. plotted in units of N /(rleQV c ) as a 

function of the photon energy ratio e/eo for the indicated values of the spatial variable y. In 
this example we have set the absorption constant (3 — 0.4 and the source location parameter 
yo = 0.9, so that the source is located near the base of the accretion column. 



This formula provides an extremely efficient alternative to numerical integration for the computation of <L 



3.5 Numerical Examples 



In this section we illustrate the computational method by examining the dependence of the Green's 
function f G (yo,y,£o,e) on the spatial location y and the energy e. We remind the reader that the solution 
for the Green's function represents the photon spectrum inside the accretion column at the specified position 
and energy, resulting from the injection of monochromatic photons with energy eo from a source located 
at yo. Hence analysis of f G allows us to explore the competing effects of Fermi energization and diffusion 
as photons travel through the column. The Green's function can be computed by combining l|49() . I|51|l . 
and (|69[l once the eigenvalues A„ have been determined using (|44|l . The eigenfunction expansion for f G 
converges fairly rapidly, and in general one obtains at least five decimal digits of accuracy if the series in 
(|49l) is terminated after the first 20 terms. 

The Green's function f G (yo,y,£o^) is plotted as a function of the energy ratio e/eo & n d the location y 
in Figure 2 for the parameter values [3 = 0.4 and yo = 0.9. In this case the first eigenvalue is given by 
Ao = 4.231 (see Fig. 1), which is equal to the high-energy slope of the Green's function in the log-log plots 
in Fig. 2. The selected value of yo corresponds to a source located near the bottom of the accretion column, 
just above the stellar surface. At the source location, y = yo = 0.9, the energy spectrum extends down to 
the injection energy, eo- However, at all other radii the spectrum displays a steep turnover above that energy 
because all of the photons have experienced Fermi energization due to collisions with the infalling electrons. 
The photons with energy e — eo at the source location have been injected so recently that they have not yet 
experienced significant energization. Note that in the far upstream region (i.e., for small values of y), the 
spectrum is greatly attenuated due to the inability of the photons to diffuse upstream through the rapidly 
infalling plasma. In this example, the average photon energy achieves its maximum value in the upstream 
region because these are the photons that have resided in the flow the longest and therefore experienced the 
most energy amplification. However, due to the attenuation mentioned above, there are not many of these 
photons. 
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Figure 3. Same as Fig. 2, except (3 — 0.4 and yo — 0.9. In this case the source is located 
in the upstream region, and the average photon energy achieves its maximum value in the 
downstream region. 

In Fig. 3 we plot the Green's function / G for the case with (3 — 4 and yo — 0.4, which yields for the first 
eigenvalue Ao = 6.325. The source is now located in the upstream region and the absorption is stronger, and 
consequently the behavior is somewhat different from that displayed in Fig. 2. In particular, the photons 
experience less overall compression in the flow and therefore the spectrum is steeper at high energies, as 
evidenced by the increase in the primary eigenvalue Ao- This is mainly due to the larger value of (3, which 
causes the photons to spend less time on average in the flow being energized by collisions with the electrons 
before they escape from the column or are "recycled" by absorption. We also note that in this case the average 
radiation energy displays its maximum value in the downstream region. This is the reverse of the behavior 
displayed in Fig. 2 because in the present situation, the source is located in the upstream region and therefore 
the photons that diffuse further upstream do not experience as much energization as those considered in 
Fig. 2. The radiation distribution in the far upstream region is greatly attenuated due to diffusion against 
the current of infalling electrons, as in Fig. 2. The analytical results for the Green's function obtained here 
provide the basis for the consideration of any source distribution since the fundamental differential equation 
Q is linear. This is further discussed in section VI. 



We can derive two interesting summation formulas for the hypergeometric eigenfunctions by using the 
transport equation 1)11(1 to study the behavior of the "energy moments," It, defined by 



The lower bound of e is chosen because / G = for e < e as explained in the discussion preceding fT3*|l . 
Note that according to 10 , the number and energy densities are given by n G = I2 and U G = I3, respectively. 
The differential equation satisfied by It is obtained by operating on 1(11(1 with J e e de, which yields 



IV. HYPERGEOMETRIC SUMMATION FORMULA 




(71) 



2/(1 - y) 



d 2 i, 



dy 2 



t 



+ 



( 



l-5y \ die 
4 J dy 
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3(3v 6(y-y )I e 3 iV e^ 2 S(y - y ) 
7v c 7-Krlvc 



(72) 



The energy moment Ig must be continuous at y = yo in order to avoid generating an infinite spatial diffusion 
flux there, and consequently we have 



A [h{y)\ 



. (73) 



y=Va 

By integrating (|72|l in a small region around y = yo, we can show that Ig also satisfies the derivative jump 
condition 



A 



dig 



dy 



3 0It(Vo) _ ZNqe^ 2 (?4) 



4y l-Kr^v c y {l-y Q ) 

y=yo 

where we have also utilized (|12ll . 

The homogeneous version of (|72|l obtained when y ^ yo is equivalent to Ijl4(l for g if we replace A with 
£ + 1. Since the energy moments Ig must satisfy the same upstream and downstream boundary conditions 
that apply to the separation eigenfunctions g n , we can therefore write the general solution for Ig as 

m = i Ceipii * +1 ' y) ' y - yo ' (75) 

[Dgip 2 {l + l,y) , y>y , 

where the constants Cg and Dg are computed by satisfying the continuity and derivative jump conditions 
given by (|75|l and Upon substitution, we obtain after some algebra 



77rv c rl 3/3 pi (* + l,i/o) <P2{* + l,Vo) -4 y W(£ + l,y ) ' 



_ I2iy 0£ ^ 2 (i-yo)-Vi(l + i,yo) f77 s 

where W(f + l,yo) is computed using [cf. Ea.l|4"5)l] 

-1/4 



^ + ^) = i rr r m q ; ) /^ f— » (78) 

4 r(a^)r(-l/4) 1 - y 



and 



9 - V33+ 16 £ 



(79) 



The energy moments Ig(y) can also be calculated by substituting for the Green's function in the 
fundamental integral (|71|l using l|49|) . Reversing the order or summation and integration yields 

oo 

Mv) = 4 +1 T, A n {X n -l-l)- 1 9n {y) , (80) 

71=0 

where g n (y) and A n are given by (|27|l and l|51|l . respectively. Note that the expression for g n (y) can be 
rewritten as 

_ (/3l(A n ,y mi n) V2(A 
9n{y) = 77 \ ! y° L ) 

where 

2/min = min(y, y ) , y max = max(y, y ) . (82) 

Eliminating Ig(y) between (|75|l and H8U|I and making use of l)5ip. I|7t)l) . (|77|l . and (|81|) . we find after some 
simplification that 

E Pi (A n , yo ) Pi(A„,y min )p2(A ni ymax) 
7 (l-yo) _1 Pi(^ + l,y mi „)p2^ + l j ymax J /oo\ 
=o p 2 (A„,y ) (A„ - i - 1) ~ 3/3 ^ (£ + 1, y ) p 2 (£ + 1, y ) - Ay W(£ + 1, yo) ' 1 ' 

where the eigenvalues X n are computed using (jUjl . Equation (JHSJ is a new hypergeometric summation 
formula that has not appeared previously in the literature. This relation holds for all real values of £. 
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V. LINEAR AND BILINEAR GENERATING FUNCTIONS 



The case with (3 — is interesting from a mathematical point of view because in this limit, the 
hypergeometric eigenfunctions reduce to Jacobi polynomials. We can therefore combine various results from 
sections III and IV to obtain two new summation formulas (i.e., linear and bilinear generating functions) 
for the Jacobi polynomials that have not appeared previously in the literature. In the limit fi — > 0, the 
eigenvalue equation (pHJl reduces to 

W(A n , yo) = ~ r( T( ;i; ^ = , (84) 

4 r(a)r(-l/4) 1 - yo 

where we have also made use of (|43|l . Roots of this expression occur where |T(a)| — > oo, which corresponds 
to 

a = -n, n = 0,1,2,... (85) 
In this situation, we can use l|19|) to demonstrate that the exact solution for the eigenvalues X n is given by 

A„ = 4n 2 + 9n + 4. (86) 
Next we note that a + b — 9/4 in general according to (|19fl . and therefore we find that 

b=l+n. (87) 

The corresponding expression for the fundamental upstream eigensolution, ipi(X n ,y), is given in this case by 
the polynomial [see Eq. (|17|) ] 

/ 9 9 
Vi{K,y) = yF l-n, - + n; y 

and the fundamental eigensolution in the downstream region, <f2(X n ,y), likewise reduces to [see Eq. (|25|) ] 

Hence the two eigensolutions <pi(X n ,y) and </?2(A n ,y) are linearly dependent functions in this case, which is 
expected since the Wronskian W(X n ,yo) = according to 184fl . This in turn reflects the fact that there is 
no derivative jump in the global separation eigenfunction g n {y) at y = yo when (3 — [see Eq. <|16|) ]. 

Due to the linear dependence of t/?i(A n , y) and <p2(A„, y), equation l|81l) for the global eigenfunction g n (y) 
now simplifies to 

9n(y) = (fx(X n ,y) , (90) 
and therefore the summation formula presented in (|83|l can be rewritten in the j3 = case as 

(An-^-l)Cn = Ay^(l~y )W(e+l,y ) ' 

where y min and y max are defined by and W(£ + 1, yo) is computed using (|7S|l . 

We are now in a position to derive an interesting summation formula for products of Jacobi polynomials. 
Using Eq. (15.4.6) from Abramowitz & Stegunpl our expression for the eigensolution <,5i(A n ,y) can be 
rewritten as 

vAK,y) = j^jt »4 5/4 ' 0> (l-2») , (92) 

pW-d-^^F^^;?;,) (93) 
represents the Jacobi polynomial, and (a) n denotes the Pochhammer symbol, defined b}0 
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where 



In the present application, with (3 — 0, we can combine 11521 . (|90|l . and (|92|) to express the quadratic 
normalization integrals, C„, as 



n 



2 ,1 



y 5/i 



2 



Pi 5 / 4 '°)(l-2 2/ ) dy, (95) 



.(9/4)„j j o 

which can be evaluated using Eq. (7.391.1) from Gradshteyn and Ryzhil^lto obtain 

n\ 



£ = 



gx-l 



2n +i ) . (96) 



.(9/4) 

Equations JZHJl, ®, (EU, and JHH) can be combined to derive a new bilinear generating function for 

the Jacobi polynomials, which can be written as 

f fa , o n Pi 5/4 ' 0) (l-2y )pl 5/4 ' 0) (l-2^) _ 16 r(3/4)r(a<) + 1, y min ) y 2 (l + 1, g/ max ) 

+ J 4n2 + 9n + 3-£ "5 r(l - o< ) yy ' 1 J 

where ai is defined by (|79|l . Note that the functions ipi(£ + l,y m in) and + 1,2/max) appearing on the 
right-hand side of l|97|) are not eigenfunctions since in general the quantity I + 1 is not equal to one of the 
eigenvalues A„. 

An interesting special case occurs in the limit yo — * 0. Making use of the relation [see Eq. (|93|l ] 

p(5/4,0) (1)= M ! » ) (9g) 

and the identity 

r (!) r (?)=iH 1/2 . <"» 

we now find that i|97|) reduces to the linear generating function 

i / 4 9 :rf+3 + T' ^ 5, " , < i - 2 »» - • < io »' 

(4n2 + 9n + 3 — l) n\ 1(1 — ai) y 

which is valid for all real values of £. Equations Q97JI and 1)100(1 are new results that are useful for the 
evaluation of infinite sums containing either products of Jacobi polynomials or single Jacobi polynomials, 
respectively. 

VI. CONCLUSION 

In this article we have employed methods of classical analysis to obtain the exact solution for the Green's 
function describing the Fermi energization of photons scattered by infalling electrons in a pulsar accretion 
column. This process is of central importance in the development of theoretical models for the production 
of the X-ray spectra observed from these objects, which are among the brightest sources in the Milky Way 
galaxy. As demonstrated in Fig. 1 and equation Q49JI. the Green's function is characterized by a power- 
law shape at high photon energies, which is typical for a Fermi process. In this scenario, photons gain 
their energy by diffusing back and forth across the shock many times. The probability of multiple shock 
crossings decreases exponentially with the number of crossings, and the mean energy of the photons increases 
exponentially with the number of crossings. This combination of factors naturally gives rise to a power-law 
energy distribution.UJHence shock energization in the pulsar accretion column provides a natural explanation 
for the spectrum of the high-energy radiation produced by X-ray pulsars. Specific examples of the Green's 
function are plotted in Figs. 2 and 3. 

Due to the linearity of the transport equation (Q, we can employ the Green's function to calculate 
the radiati on s pectrum inside the accretion column resulting from an arbitrary source spectrum using the 
convolutiorll^l 

/U/o,y,e)=/ j{to)— £ de , (101) 

Jo iv 

where j(eo)cfeo represents the number of photons injected into the accretion column between at location 
yo with energy between eo and eo + deo- The source distribution of greatest astrophysical interest is the 
"thermal mound" source located near the base of the accretion column, where the gas has decelerated 
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almost to rest and is therefore extremely dense. This hot plasma is in full thermodynamic equilibrium, 
and consequently it radiates a blackbody spectrumP The absorption parameter (3 has been included in the 
transport equation in order to account for the fact that a blackbody acts as both a source and a sink of 
radiationPl The fundamental results for the Green's function obtained in the present article will be used to 
study the reprocessing of the blackbody radiation emitted from the thermal mound in a subsequent paper. 

In addition to the analytical results for the Green's function, we have also obtained an interesting formula 
for the evaluation of an infinite series involving products of the orthogonal hypergeometric eigenfunctions 
[see Eq. (|83H ]. This derivation was based on the simultaneous calculation of the energy moments Ig(y) 
using either an expression based on term- by-term integration of the Green's function expansion (|49|l , or an 
independent solution developed via direct integration of the fundamental transport equation In the 
special case (3 — > 0, which corresponds physically to the neglect of absorption at the source location, our 
general formula for the hypergemetric summation reduces to a bilinear generating function for the Jacobi 
polynomials given by (|97|) . This relation in turn simplifies to yield a linear generating function for the Jacobi 
polynomials in the limit ya — > 0, which corresponds physically to a source located in the far upstream region 
[see Eq. (ITUrtfl ]. 

The results derived is this article for the linear and bilinear generating functions of Jacob i polynomi als 
are related to v arious similar expressions obtained previously by Chen and SrivastavaplEl Srivastavap^l 
Rangarajanp^l anc j Pittaluga, Sacripante, and Srivastava!^ However, our results are not identical to any 
of their formulas and therefore they represent an interesting new family of relations. Although the linear 
and bilinear generating functions developed here relate specifically to the properties of the polynomials 
p^(5/4, o)^ _ we ex p ec t that some level of generalization may be possible. We plan to pursue this 
question in future work. 
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